
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%                 TABLE OS.3:               %%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')





pretemp = pre ;
pre     = pre*36 ;


format shortg
disp('*****************************************************************************************')
disp('*****************************************************************************************')
disp('****************    INITIAL PROFICIENCY PRODUCTION FUNCTION (RF):   *********************')
disp('*****************************************************************************************')
disp('*****************************************************************************************')

logPreScore = log(pre) ;  
labels_LR1  = {'Constant';'log(thetaE)';'log(thetaL)'} ;
ExpVars_LR1 = [ones(size(logPreScore)), ...
               log(Y_EimputeShrunk), ...
               log(Y_LimputeShrunk)] ; 
[alpha_LR1,se_LR1,tstat_LR1,pval_LR1,CI_LR1,out_LR1] = regressHet(logPreScore,ExpVars_LR1,'FGLSrobust',labels_LR1,[],[],0,[],[],[]); 
alphaE_LR1  = alpha_LR1(2)*ones(length(logPreScore),1) ;  
alphaL_LR1  = alpha_LR1(3)*ones(length(logPreScore),1) ; 
TFP_LR1     = exp( alpha_LR1(1)*ones(length(logPreScore),1) ) ; 
StDevlogPre        = std(logPreScore) ;
StDevEffectATE     = nan(length(alpha_LR1)+2,1) ;
StDevEffectATE(2)  = std(log(Y_EimputeShrunk))*alpha_LR1(2)/StDevlogPre ;  %%%%This is a one standard deviation increase in log(thetaE)
StDevEffectATE(3)  = std(log(Y_LimputeShrunk))*alpha_LR1(3)/StDevlogPre ;  %%%%This is a one standard deviation increase in log(thetaL)
disp('-----------------------------------------------------------------------------------')
disp('-------------------  INITIAL PROFICIENCY PRODUCTION: MODEL 1  -------------------')
disp('-----------------------------------------------------------------------------------')
out_LR1.results.StDevEffectATE = StDevEffectATE;
WhiteTest  = out_LR1.White %#ok
Method     = out_LR1.method %#ok
disp(out_LR1.results(:,:)) 

%%

labels_LR2  = {'Constant';'log(thetaE)';'log(thetaL)';'D2';'D3';'D2*log(thetaE)';'D3*log(thetaE)';'D2*log(thetaL)';'D3*log(thetaL)'} ;
ExpVars_LR2 = [ones(size(logPreScore)), ...1
               log(Y_EimputeShrunk), ...2
               log(Y_LimputeShrunk), ...3
               D2, ...4
               D3, ...5
               log(Y_EimputeShrunk).*D2, ...6
               log(Y_EimputeShrunk).*D3, ...7
               log(Y_LimputeShrunk).*D2, ...8
               log(Y_LimputeShrunk).*D3] ; %9
[alpha_LR2,se_LR2,tstat_LR2,pval_LR2,CI_LR2,out_LR2] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,[],[],0,[],[],[]); 
alphaE_LR2  = alpha_LR2(2)*ones(length(logPreScore),1) + D2*alpha_LR2(6) + D3*alpha_LR2(7) ;  
alphaL_LR2  = alpha_LR2(3)*ones(length(logPreScore),1) + D2*alpha_LR2(8) + D3*alpha_LR2(9) ;
TFP_LR2     = exp( alpha_LR2(1)*ones(length(logPreScore),1) + D2*alpha_LR2(4) + D3*alpha_LR2(5) ) ;  
StDevEffectATE = nan(length(alpha_LR2)+2,1) ;
StDevEffectATE(1) = std( log(TFP_LR2) )/StDevlogPre ; %%%%This represents a one standard deviation change in the log(TFP) factor, which varies across individuals
StDevEffectATE(2) = mean( std(log(Y_EimputeShrunk))*alphaE_LR2 )/StDevlogPre  ;  %%%%This a one standard deviation increase in log(thetaE)
StDevEffectATE(3) = mean( std(log(Y_LimputeShrunk))*alphaL_LR2 )/StDevlogPre ;  %%%%This a one standard deviation increase in log(thetaL)
StDevEffectATE(4) = mean( ( alpha_LR2(4) + log(Y_EimputeShrunk)*alpha_LR2(6) + log(Y_LimputeShrunk)*alpha_LR2(8) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 2
StDevEffectATE(5) = mean( ( alpha_LR2(5) + log(Y_EimputeShrunk)*alpha_LR2(7) + log(Y_LimputeShrunk)*alpha_LR2(9) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 3
StDevEffectAToT = nan(length(alpha_LR2)+2,1) ;
StDevEffectAToT(4) = mean( ( alpha_LR2(4) + log(Y_EimputeShrunk(D2==1))*alpha_LR2(6) + log(Y_LimputeShrunk(D2==1))*alpha_LR2(8) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 2
StDevEffectAToT(5) = mean( ( alpha_LR2(5) + log(Y_EimputeShrunk(D3==1))*alpha_LR2(7) + log(Y_LimputeShrunk(D3==1))*alpha_LR2(9) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 3
disp('-----------------------------------------------------------------------------------')
disp('-------------------  INITIAL PROFICIENCY PRODUCTION: MODEL 2  -------------------')
disp('-----------------------------------------------------------------------------------')
out_LR2.results.StDevEffectATE = StDevEffectATE;
out_LR2.results.StDevEffectAToT = StDevEffectAToT;
WhiteTest  = out_LR2.White %#ok
Method     = out_LR2.method %#ok
disp(out_LR2.results(:,:)) 
components = {'log(TFP)';'alphaE';'alphaL'}; Mean = [mean(log(TFP_LR2)); mean(alphaE_LR2); mean(alphaL_LR2)]; StDev = [std(log(TFP_LR2)); std(alphaE_LR2); std(alphaL_LR2)] ; 
MeanStDevAcrossi = table(components,Mean,StDev)  %#ok
%%%%Test for joint significance of TFP:
R = zeros(3,length(ExpVars_LR2(1,:))) ; R(1,1)=1; R(2,4)=1; R(3,5)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
TFPtest          = temp.Wald %#ok
%%%%Test for joint significance of alphaE terms:
R = zeros(3,length(ExpVars_LR2(1,:))) ; R(1,2)=1; R(2,6)=1; R(3,7)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
logEtest   = temp.Wald %#ok
%%%%Test for joint significance of alphaL terms:
R = zeros(3,length(ExpVars_LR2(1,:))) ; R(1,3)=1; R(2,8)=1; R(3,9)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
logLtest   = temp.Wald %#ok
%%%%Test for joint significance of D2 terms:
R = zeros(3,length(ExpVars_LR2(1,:))) ; R(1,4)=1; R(2,6)=1; R(3,8)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
D2test   = temp.Wald %#ok
%%%%Test for joint significance of D2 SLOPE terms:
R = zeros(2,length(ExpVars_LR2(1,:))) ; R(1,6)=1; R(2,8)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
D2Slopetest   = temp.Wald %#ok
%%%%Test for joint significance of D3 terms:
R = zeros(3,length(ExpVars_LR2(1,:))) ; R(1,5)=1; R(2,7)=1; R(3,9)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
D3test   = temp.Wald %#ok
%%%%Test for joint significance of D3 SLOPE terms:
R = zeros(2,length(ExpVars_LR2(1,:))) ; R(1,7)=1; R(2,9)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
D3Slopetest   = temp.Wald %#ok
%%%%Test for joint significance of ALL D2&D3 terms:
R = zeros(6,length(ExpVars_LR2(1,:))) ; for ii=1:6; R(ii,3+ii)=1; end; r = zeros(6,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
D23test   = temp.Wald %#ok
%%%%Test for joint significance of D2&D3 SLOPE terms:
R = zeros(4,length(ExpVars_LR2(1,:))) ; R(1,6)=1; R(2,7)=1; R(3,8)=1; R(4,9)=1; r = zeros(4,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
D23Slopetest   = temp.Wald %#ok
%%%%Test for joint significance of D2&D3 TFP terms:
R = zeros(2,length(ExpVars_LR2(1,:))) ; R(1,4)=1; R(2,5)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR2,'FGLSrobust',labels_LR2,R,r,0,[],[],[]); 
D23TFPtest   = temp.Wald %#ok



%%


labels_LR3  = {'Constant';'log(thetaE)';'log(thetaL)';'D2';'D3';'D2*log(thetaE)';'D3*log(thetaE)';'D2*log(thetaL)';'D3*log(thetaL)';...
               'grade5';'grade5*log(thetaE)';'grade5*log(thetaL)'; 'fem'; 'fem*log(thetaE)'; 'fem*log(thetaL)';...
               'black';'hispanic';'log(thetaE)*blk';'log(thetaE)*hsp';'log(thetaL)*blk';'log(thetaL)*hsp'} ;
ExpVars_LR3 = [ones(size(logPreScore)), ...1
               log(Y_EimputeShrunk), ...2
               log(Y_LimputeShrunk), ...3
               D2, ...4
               D3, ...5
               log(Y_EimputeShrunk).*D2, ...6
               log(Y_EimputeShrunk).*D3, ...7
               log(Y_LimputeShrunk).*D2, ...8
               log(Y_LimputeShrunk).*D3, ...9
               grade5, ...10
               grade5.*log(Y_EimputeShrunk), ...11
               grade5.*log(Y_LimputeShrunk), ...12
               fem, ...13
               fem.*log(Y_EimputeShrunk), ...14
               fem.*log(Y_LimputeShrunk), ...15
               blk, ... 16
               hsp, ... 17
               log(Y_EimputeShrunk).*blk, ... 18
               log(Y_EimputeShrunk).*hsp, ... 19
               log(Y_LimputeShrunk).*blk, ... 20
               log(Y_LimputeShrunk).*hsp] ; %21
[alpha_LR3,se_LR3,tstat_LR3,pval_LR3,CI_LR3,out_LR3] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,[],[],0,[],[],[]); 
alphaE_LR3  = alpha_LR3(2)*ones(length(logPreScore),1) + D2*alpha_LR3(6) + D3*alpha_LR3(7) + grade5*alpha_LR3(11) + fem*alpha_LR3(14) + blk*alpha_LR3(18) + hsp*alpha_LR3(19) ;  
alphaL_LR3  = alpha_LR3(3)*ones(length(logPreScore),1) + D2*alpha_LR3(8) + D3*alpha_LR3(9) + grade5*alpha_LR3(12) + fem*alpha_LR3(15) + blk*alpha_LR3(20) + hsp*alpha_LR3(21) ; 
TFP_LR3     = exp( alpha_LR3(1)*ones(length(logPreScore),1) + D2*alpha_LR3(4) + D3*alpha_LR3(5) + grade5*alpha_LR3(10) + fem*alpha_LR3(13) + blk*alpha_LR3(16) + hsp*alpha_LR3(17) ) ;  
StDevEffectATE = nan(length(alpha_LR3)+2,1) ;
StDevEffectATE(1)  = std( log(TFP_LR3) )/StDevlogPre ; %%%%This represents a one standard deviation change in the log(TFP) factor, which varies across individuals
StDevEffectATE(2)  = mean( std(log(Y_EimputeShrunk))*alphaE_LR3 )/StDevlogPre ;  %%%%This a one standard deviation increase in log(thetaE)
StDevEffectATE(3)  = mean( std(log(Y_LimputeShrunk))*alphaL_LR3 )/StDevlogPre ;  %%%%This a one standard deviation increase in log(thetaL)
StDevEffectATE(4)  = mean( ( alpha_LR3(4) + log(Y_EimputeShrunk)*alpha_LR3(6) + log(Y_LimputeShrunk)*alpha_LR3(8) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 2
StDevEffectATE(5)  = mean( ( alpha_LR3(5) + log(Y_EimputeShrunk)*alpha_LR3(7) + log(Y_LimputeShrunk)*alpha_LR3(9) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 3
StDevEffectATE(10) = mean( ( alpha_LR3(10) + log(Y_EimputeShrunk)*alpha_LR3(11) + log(Y_LimputeShrunk)*alpha_LR3(12) )/StDevlogPre ) ;  %%%%This is a change from grade 6 to grade 5
StDevEffectATE(13) = mean( ( alpha_LR3(13) + log(Y_EimputeShrunk)*alpha_LR3(14) + log(Y_LimputeShrunk)*alpha_LR3(15) )/StDevlogPre ) ;  %%%%This is a change from male to female
StDevEffectATE(16) = mean( ( alpha_LR3(16) + log(Y_EimputeShrunk)*alpha_LR3(18) + log(Y_LimputeShrunk)*alpha_LR3(20) )/StDevlogPre ) ;  %%%%This is a change from WAO to BLACK
StDevEffectATE(17) = mean( ( alpha_LR3(17) + log(Y_EimputeShrunk)*alpha_LR3(19) + log(Y_LimputeShrunk)*alpha_LR3(21) )/StDevlogPre ) ;  %%%%This is a change from WAO to HISPANIC
StDevEffectAToT = nan(length(alpha_LR3)+2,1) ;
StDevEffectAToT(4)  = mean( ( alpha_LR3(4) + log(Y_EimputeShrunk(D2==1))*alpha_LR3(6) + log(Y_LimputeShrunk(D2==1))*alpha_LR3(8) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 2
StDevEffectAToT(5)  = mean( ( alpha_LR3(5) + log(Y_EimputeShrunk(D3==1))*alpha_LR3(7) + log(Y_LimputeShrunk(D3==1))*alpha_LR3(9) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 3
StDevEffectAToT(10) = mean( ( alpha_LR3(10) + log(Y_EimputeShrunk(grade5==1))*alpha_LR3(11) + log(Y_LimputeShrunk(grade5==1))*alpha_LR3(12) )/StDevlogPre ) ;  %%%%This is a change from grade 6 to grade 5
StDevEffectAToT(13) = mean( ( alpha_LR3(13) + log(Y_EimputeShrunk(fem==1))*alpha_LR3(14) + log(Y_LimputeShrunk(fem==1))*alpha_LR3(15) )/StDevlogPre ) ;  %%%%This is a change from male to female
StDevEffectAToT(16) = mean( ( alpha_LR3(16) + log(Y_EimputeShrunk(blk==1))*alpha_LR3(18) + log(Y_LimputeShrunk(blk==1))*alpha_LR3(20) )/StDevlogPre ) ;  %%%%This is a change from WAO to BLACK
StDevEffectAToT(17) = mean( ( alpha_LR3(17) + log(Y_EimputeShrunk(hsp==1))*alpha_LR3(19) + log(Y_LimputeShrunk(hsp==1))*alpha_LR3(21) )/StDevlogPre ) ;  %%%%This is a change from WAO to HISPANIC
disp('-----------------------------------------------------------------------------------')
disp('-------------------  INITIAL PROFICIENCY PRODUCTION: MODEL 3  -------------------')
disp('-----------------------------------------------------------------------------------')
out_LR3.results.StDevEffectATE = StDevEffectATE;
out_LR3.results.StDevEffectAToT = StDevEffectAToT;
WhiteTest  = out_LR3.White %#ok
Method     = out_LR3.method %#ok
disp(out_LR3.results(:,:)) 
components = {'log(TFP)';'alphaE';'alphaL'}; Mean = [mean(log(TFP_LR3)); mean(alphaE_LR3); mean(alphaL_LR3)]; StDev = [std(log(TFP_LR3)); std(alphaE_LR3); std(alphaL_LR3)] ; 
MeanStDevAcrossi = table(components,Mean,StDev) %#ok
%%%%Test for joint significance of TFP:
R = zeros(7,length(ExpVars_LR3(1,:))) ; R(1,1)=1; R(2,4)=1; R(3,5)=1; R(4,10)=1; R(5,13)=1; R(6,16)=1; R(7,17)=1; r = zeros(7,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
TFPtest          = temp.Wald %#ok
%%%%Test for joint significance of alphaE terms:
R = zeros(7,length(ExpVars_LR3(1,:))) ; R(1,2)=1; R(2,6)=1; R(3,7)=1; R(4,11)=1; R(5,14)=1; R(6,18)=1; R(7,19)=1; r = zeros(7,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
logEtest   = temp.Wald %#ok
%%%%Test for joint significance of alphaL terms:
R = zeros(7,length(ExpVars_LR3(1,:))) ; R(1,3)=1; R(2,8)=1; R(3,9)=1; R(4,12)=1; R(5,15)=1; R(6,20)=1; R(7,21)=1; r = zeros(7,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
logLtest   = temp.Wald %#ok
%%%%Test for joint significance of D2 terms:
R = zeros(3,length(ExpVars_LR3(1,:))) ; R(1,4)=1; R(2,6)=1; R(3,8)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
D2test   = temp.Wald %#ok
%%%%Test for joint significance of D2 SLOPE terms:
R = zeros(2,length(ExpVars_LR3(1,:))) ; R(1,6)=1; R(2,8)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
D2Slopetest   = temp.Wald %#ok
%%%%Test for joint significance of D3 terms:
R = zeros(3,length(ExpVars_LR3(1,:))) ; R(1,5)=1; R(2,7)=1; R(3,9)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
D3test   = temp.Wald %#ok
%%%%Test for joint significance of D3 SLOPE terms:
R = zeros(2,length(ExpVars_LR3(1,:))) ; R(1,7)=1; R(2,9)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
D3Slopetest   = temp.Wald %#ok
%%%%Test for joint significance of ALL D2&D3 terms:
R = zeros(6,length(ExpVars_LR3(1,:))) ; for ii=1:6; R(ii,3+ii)=1; end; r = zeros(6,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
D23test   = temp.Wald %#ok
%%%%Test for joint significance of D2&D3 SLOPE terms:
R = zeros(4,length(ExpVars_LR3(1,:))) ; R(1,6)=1; R(2,7)=1; R(3,8)=1; R(4,9)=1; r = zeros(4,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
D23Slopetest   = temp.Wald %#ok
%%%%Test for joint significance of D2&D3 TFP terms:
R = zeros(2,length(ExpVars_LR3(1,:))) ; R(1,4)=1; R(2,5)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
D23TFPtest   = temp.Wald %#ok
%%%%Test for joint significance of grade5 terms:
R = zeros(3,length(ExpVars_LR3(1,:))) ; R(1,10)=1; R(2,11)=1; R(3,12)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
grade5test   = temp.Wald %#ok
%%%%Test for joint significance of gender terms:
R = zeros(3,length(ExpVars_LR3(1,:))) ; R(1,13)=1; R(2,14)=1; R(3,15)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
femtest   = temp.Wald %#ok
%%%%Test for joint significance of blk terms:
R = zeros(3,length(ExpVars_LR3(1,:))) ; R(1,16)=1; R(2,18)=1; R(3,20)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
blktest   = temp.Wald %#ok
%%%%Test for joint significance of hsp terms:
R = zeros(3,length(ExpVars_LR3(1,:))) ; R(1,17)=1; R(2,19)=1; R(3,21)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR3,'FGLSrobust',labels_LR3,R,r,0,[],[],[]); 
hsptest   = temp.Wald %#ok


% NHelpers    = AchievementData.helpadult + AchievementData.helppeer ; NHelpers = NHelpers(useindx) ;
% NHelppeer   = AchievementData.helppeer ;  NHelppeer  = NHelppeer(useindx) ;
% NHelpadult  = AchievementData.helpadult ; NHelpadult = NHelpadult(useindx) ;
% Helperuseindx = find(~isnan(NHelpers)) ;
NHelpers    = X_E(:,4) + X_E(:,5) ; NHelpers = NHelpers(useindx) ;
NHelppeer   = X_E(:,5) ; NHelppeer  = NHelppeer(useindx) ;
NHelpadult  = X_E(:,4) ; NHelpadult = NHelpadult(useindx) ;
NoHmIntrnt  = X_E(:,12) ; 
MoblFrac    = X_E(:,13) ; 
TabltFrac   = X_E(:,14) ;
%%

labels_LR4  = {'Constant';'log(thetaE)';'log(thetaL)';'D2';'D3';'D2*log(thetaE)';'D3*log(thetaE)';'D2*log(thetaL)';'D3*log(thetaL)';...
               'grade5';'grade5*log(thetaE)';'grade5*log(thetaL)'; 'fem'; 'fem*log(thetaE)'; 'fem*log(thetaL)';...
               'black';'hispanic';'log(thetaE)*blk';'log(thetaE)*hsp';'log(thetaL)*blk';'log(thetaL)*hsp';...
               'SSEInc';'SSEIns';'log(thetaE)*SSEInc';'log(thetaE)*SSEIns';'log(thetaL)*SSEInc';'log(thetaL)*SSEIns';...
               'NHelpers';'log(thetaE)*NHelpers';'log(thetaL)*NHelpers';...
               'NoHmIntrnt';'log(thetaE)*NoHmIntrnt';'log(thetaL)*NoHmIntrnt';...
               'MoblFrac';'log(thetaE)*MoblFrac';'log(thetaL)*MoblFrac';...
               'TabltFrac';'log(thetaE)*TabltFrac';'log(thetaL)*TabltFrac'} ;
ExpVars_LR4 = [ones(size(logPreScore)), ...1
               log(Y_EimputeShrunk), ...2
               log(Y_LimputeShrunk), ...3
               D2, ...4
               D3, ...5
               log(Y_EimputeShrunk).*D2, ...6
               log(Y_EimputeShrunk).*D3, ...7
               log(Y_LimputeShrunk).*D2, ...8
               log(Y_LimputeShrunk).*D3, ...9
               grade5, ...10
               grade5.*log(Y_EimputeShrunk), ...11
               grade5.*log(Y_LimputeShrunk), ... 12
               fem, ...13
               fem.*log(Y_EimputeShrunk), ...14
               fem.*log(Y_LimputeShrunk), ...15
               blk, ... 16
               hsp, ... 17
               log(Y_EimputeShrunk).*blk, ... 18
               log(Y_EimputeShrunk).*hsp, ... 19
               log(Y_LimputeShrunk).*blk, ... 20
               log(Y_LimputeShrunk).*hsp, ...21
               SSEInc, ... 22
               SSEIns, ... 23
               log(Y_EimputeShrunk).*SSEInc, ... 24
               log(Y_EimputeShrunk).*SSEIns, ... 25
               log(Y_LimputeShrunk).*SSEInc, ... 26
               log(Y_LimputeShrunk).*SSEIns, ... 27 
               NHelpers, ... 28 
               log(Y_EimputeShrunk).*NHelpers, ... 29
               log(Y_LimputeShrunk).*NHelpers, ... 30
               NoHmIntrnt, ... 31
               log(Y_EimputeShrunk).*NoHmIntrnt, ... 32
               log(Y_LimputeShrunk).*NoHmIntrnt, ... 33
               MoblFrac, ... 34
               log(Y_EimputeShrunk).*MoblFrac, ... 35
               log(Y_LimputeShrunk).*MoblFrac, ... 36
               TabltFrac, ... 37
               log(Y_EimputeShrunk).*TabltFrac, ... 38
               log(Y_LimputeShrunk).*TabltFrac] ; % 39
%%%%Test for joint significance of all the SSEInc & SSEIns terms
R = zeros(9,length(ExpVars_LR4(1,:))) ; for ii=1:9; R(ii,21+ii) = 1 ; end; r = zeros(9,1) ;
[alpha_LR4,se_LR4,tstat_LR4,pval_LR4,CI_LR4,out_LR4] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]);  
alphaE_LR4  = alpha_LR4(2)*ones(length(logPreScore),1) + D2*alpha_LR4(6) + D3*alpha_LR4(7) + grade5*alpha_LR4(11) ...
              + fem*alpha_LR4(14) + blk*alpha_LR4(18) + hsp*alpha_LR4(19) + SSEInc*alpha_LR4(24) + SSEIns*alpha_LR4(25) ...
              + NHelpers*alpha_LR4(29) + NoHmIntrnt*alpha_LR4(32) + MoblFrac*alpha_LR4(35) + TabltFrac*alpha_LR4(38) ;  
alphaL_LR4  = alpha_LR4(3)*ones(length(logPreScore),1) + D2*alpha_LR4(8) + D3*alpha_LR4(9) + grade5*alpha_LR4(12) ...
              + fem*alpha_LR4(15) + blk*alpha_LR4(20) + hsp*alpha_LR4(21) + SSEInc*alpha_LR4(26) + SSEIns*alpha_LR4(27) ...
              + NHelpers*alpha_LR4(30) + NoHmIntrnt*alpha_LR4(33) + MoblFrac*alpha_LR4(36) + TabltFrac*alpha_LR4(39) ;  
TFP_LR4     = exp( alpha_LR4(1)*ones(length(logPreScore),1) + D2*alpha_LR4(4) + D3*alpha_LR4(5) + grade5*alpha_LR4(10) ...
              + fem*alpha_LR4(13) + blk*alpha_LR4(16) + hsp*alpha_LR4(17) + SSEInc*alpha_LR4(22) + SSEIns*alpha_LR4(23) ...
              + NHelpers*alpha_LR4(28)  + NoHmIntrnt*alpha_LR4(31) + MoblFrac*alpha_LR4(34) + TabltFrac*alpha_LR4(37) ) ; 
StDevEffectATE = nan(length(alpha_LR4)+2,1) ;  
StDevEffectATE(1)  = std( log(TFP_LR4) )/StDevlogPre ; %%%%This represents a one standard deviation change in the log(TFP) factor, which varies across individuals
StDevEffectATE(2)  = mean( std(log(Y_EimputeShrunk))*alphaE_LR4 )/StDevlogPre ;  %%%%This is a one standard deviation increase in log(thetaE)
StDevEffectATE(3)  = mean( std(log(Y_LimputeShrunk))*alphaL_LR4 )/StDevlogPre ;  %%%%This is a one standard deviation increase in log(thetaL)
StDevEffectATE(4)  = mean( ( alpha_LR4(4) + log(Y_EimputeShrunk)*alpha_LR4(6) + log(Y_LimputeShrunk)*alpha_LR4(8) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 2
StDevEffectATE(5)  = mean( ( alpha_LR4(5) + log(Y_EimputeShrunk)*alpha_LR4(7) + log(Y_LimputeShrunk)*alpha_LR4(9) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 3
StDevEffectATE(10) = mean( ( alpha_LR4(10) + log(Y_EimputeShrunk)*alpha_LR4(11) + log(Y_LimputeShrunk)*alpha_LR4(12) )/StDevlogPre ) ;  %%%%This is a change from grade 6 to grade 5
StDevEffectATE(13) = mean( ( alpha_LR4(13) + log(Y_EimputeShrunk)*alpha_LR4(14) + log(Y_LimputeShrunk)*alpha_LR4(15) )/StDevlogPre ) ;  %%%%This is a change from male to female
StDevEffectATE(16) = mean( ( alpha_LR4(16) + log(Y_EimputeShrunk)*alpha_LR4(18) + log(Y_LimputeShrunk)*alpha_LR4(20) )/StDevlogPre ) ;  %%%%This is a change from WAO to BLACK
StDevEffectATE(17) = mean( ( alpha_LR4(17) + log(Y_EimputeShrunk)*alpha_LR4(19) + log(Y_LimputeShrunk)*alpha_LR4(21) )/StDevlogPre ) ;  %%%%This is a change from WAO to HISPANIC
StDevEffectATE(22) = mean( std(SSEInc)*( alpha_LR4(22) + log(Y_EimputeShrunk)*alpha_LR4(24) + log(Y_LimputeShrunk)*alpha_LR4(26) )/StDevlogPre ) ;  %%%%This is a standard deviation increase in SSEInc
StDevEffectATE(23) = mean( std(SSEIns)*( alpha_LR4(23) + log(Y_EimputeShrunk)*alpha_LR4(25) + log(Y_LimputeShrunk)*alpha_LR4(27) )/StDevlogPre ) ;  %%%%This is a standard deviation increase in SSEIns
StDevEffectATE(28) = mean( std(NHelpers)*( alpha_LR4(28) + log(Y_EimputeShrunk)*alpha_LR4(29) + log(Y_LimputeShrunk)*alpha_LR4(30) )/StDevlogPre ) ;  %%%%This is a standard deviation increase in SSEInc
StDevEffectAToT = nan(length(alpha_LR4)+2,1) ;
StDevEffectAToT(4)  = mean( ( alpha_LR4(4) + log(Y_EimputeShrunk(D2==1))*alpha_LR4(6) + log(Y_LimputeShrunk(D2==1))*alpha_LR4(8) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 2
StDevEffectAToT(5)  = mean( ( alpha_LR4(5) + log(Y_EimputeShrunk(D3==1))*alpha_LR4(7) + log(Y_LimputeShrunk(D3==1))*alpha_LR4(9) )/StDevlogPre ) ;  %%%%This is a change from District 1 to District 3
StDevEffectAToT(10) = mean( ( alpha_LR4(10) + log(Y_EimputeShrunk(grade5==1))*alpha_LR4(11) + log(Y_LimputeShrunk(grade5==1))*alpha_LR4(12) )/StDevlogPre ) ;  %%%%This is a change from grade 6 to grade 5
StDevEffectAToT(13) = mean( ( alpha_LR4(13) + log(Y_EimputeShrunk(fem==1))*alpha_LR4(14) + log(Y_LimputeShrunk(fem==1))*alpha_LR4(15) )/StDevlogPre ) ;  %%%%This is a change from male to female
StDevEffectAToT(16) = mean( ( alpha_LR4(16) + log(Y_EimputeShrunk(blk==1))*alpha_LR4(18) + log(Y_LimputeShrunk(blk==1))*alpha_LR4(20) )/StDevlogPre ) ;  %%%%This is a change from WAO to BLACK
StDevEffectAToT(17) = mean( ( alpha_LR4(17) + log(Y_EimputeShrunk(hsp==1))*alpha_LR4(19) + log(Y_LimputeShrunk(hsp==1))*alpha_LR4(21) )/StDevlogPre ) ;  %%%%This is a change from WAO to HISPANIC
disp('-----------------------------------------------------------------------------------')
disp('-------------------  INITIAL PROFICIENCY PRODUCTION: MODEL 4  -------------------')
disp('-----------------------------------------------------------------------------------')
out_LR4.results.StDevEffectATE = StDevEffectATE;
out_LR4.results.StDevEffectAToT = StDevEffectAToT;
WhiteTest  = out_LR4.White %#ok
Method     = out_LR4.method %#ok
disp(out_LR4.results(:,:)) 
components = {'log(TFP)';'alphaE';'alphaL'}; Mean = [mean(log(TFP_LR4)); mean(alphaE_LR4); mean(alphaL_LR4)]; StDev = [std(log(TFP_LR4)); std(alphaE_LR4); std(alphaL_LR4)] ; 
MeanStDevAcrossi = table(components,Mean,StDev) %#ok
%%%%Test for joint significance of TFP:
R = zeros(9,length(ExpVars_LR4(1,:))) ; R(1,1)=1; R(2,4)=1; R(3,5)=1; R(4,10)=1; R(5,13)=1; R(6,16)=1; R(7,17)=1; R(8,22)=1; R(9,23)=1; r = zeros(9,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
TFPtest          = temp.Wald %#ok
%%%%Test for joint significance of alphaE terms:
R = zeros(9,length(ExpVars_LR4(1,:))) ; R(1,2)=1; R(2,6)=1; R(3,7)=1; R(4,11)=1; R(5,14)=1; R(6,18)=1; R(7,19)=1; R(8,24)=1; R(9,25)=1; r = zeros(9,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
logEtest   = temp.Wald %#ok
%%%%Test for joint significance of alphaL terms:
R = zeros(9,length(ExpVars_LR4(1,:))) ; R(1,3)=1; R(2,8)=1; R(3,9)=1; R(4,12)=1; R(5,15)=1; R(6,20)=1; R(7,21)=1; R(8,26)=1; R(9,27)=1; r = zeros(9,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
logLtest   = temp.Wald %#ok
%%%%Test for joint significance of D2 terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,4)=1; R(2,6)=1; R(3,8)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
D2test   = temp.Wald %#ok
%%%%Test for joint significance of D2 SLOPE terms:
R = zeros(2,length(ExpVars_LR4(1,:))) ; R(1,6)=1; R(2,8)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp]  = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
D2Slopetest = temp.Wald %#ok
%%%%Test for joint significance of D3 terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,5)=1; R(2,7)=1; R(3,9)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
D3test   = temp.Wald %#ok
%%%%Test for joint significance of D3 SLOPE terms:
R = zeros(2,length(ExpVars_LR4(1,:))) ; R(1,7)=1; R(2,9)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp]  = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
D3Slopetest = temp.Wald %#ok
%%%%Test for joint significance of ALL D2&D3 terms:
R = zeros(6,length(ExpVars_LR4(1,:))) ; for ii=1:6; R(ii,3+ii)=1; end; r = zeros(6,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
D23test   = temp.Wald %#ok
%%%%Test for joint significance of D2&D3 SLOPE terms:
R = zeros(4,length(ExpVars_LR4(1,:))) ; R(1,6)=1; R(2,7)=1; R(3,8)=1; R(4,9)=1; r = zeros(4,1) ; 
[~,~,~,~,~,temp]  = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
D23Slopetest = temp.Wald %#ok
%%%%Test for joint significance of D2&D3 TFP terms:
R = zeros(2,length(ExpVars_LR4(1,:))) ; R(1,4)=1; R(2,5)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp]  = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
D23TFPtest = temp.Wald %#ok
%%%%Test for joint significance of grade5 terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,10)=1; R(2,11)=1; R(3,12)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
grade5test   = temp.Wald %#ok
%%%%Test for joint significance of grade5 interaction terms:
R = zeros(2,length(ExpVars_LR4(1,:))) ; R(1,11)=1; R(2,12)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
grade5Interactiontest   = temp.Wald %#ok
%%%%Test for joint significance of gender terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,13)=1; R(2,14)=1; R(3,15)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
femtest   = temp.Wald %#ok
%%%%Test for joint significance of gender interaction terms:
R = zeros(2,length(ExpVars_LR4(1,:))) ; R(1,14)=1; R(2,15)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp]   = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
femInteractiontest = temp.Wald %#ok
%%%%Test for joint significance of blk terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,16)=1; R(2,18)=1; R(3,20)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
blktest   = temp.Wald %#ok
%%%%Test for joint significance of blk interaction terms:
R = zeros(2,length(ExpVars_LR4(1,:))) ; R(1,18)=1; R(2,20)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
blkInteractiontest   = temp.Wald %#ok
%%%%Test for joint significance of hsp terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,17)=1; R(2,19)=1; R(3,21)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
hsptest   = temp.Wald %#ok
%%%%Test for joint significance of hsp interaction terms:
R = zeros(2,length(ExpVars_LR4(1,:))) ; R(1,19)=1; R(2,21)=1; r = zeros(2,1) ; 
[~,~,~,~,~,temp]   = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
hspInteractiontest = temp.Wald %#ok
%%%%Test for joint significance of SSEInc terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,22)=1; R(2,24)=1; R(3,26)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
SSEInctest   = temp.Wald %#ok
%%%%Test for joint significance of SSEIns terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,23)=1; R(2,25)=1; R(3,27)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
SSEInstest   = temp.Wald %#ok
%%%%Test for joint significance of SSEIns terms:
R = zeros(6,length(ExpVars_LR4(1,:))) ; for ii=1:6; R(ii,21+ii)=1; end; r = zeros(6,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
SSEtest   = temp.Wald %#ok
%%%%Test for joint significance of NHelppeer terms:
R = zeros(3,length(ExpVars_LR4(1,:))) ; R(1,28)=1; R(2,29)=1; R(3,30)=1; r = zeros(3,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
NHelperstest   = temp.Wald %#ok

%%%%Test for joint significance of Home Connectivity terms:
R = zeros(9,length(ExpVars_LR4(1,:))) ; for ii=1:9; R(ii,30+ii)=1; end; r = zeros(9,1) ; 
[~,~,~,~,~,temp] = regressHet(logPreScore,ExpVars_LR4,'FGLSrobust',labels_LR4,R,r,0,[],[],[]); 
HomeConnectivitytest   = temp.Wald %#ok
%%


%%
figure
% alphaE_LR4  = alpha_LR4(2)*ones(length(logPreScore),1) + D2*alpha_LR4(6) + D3*alpha_LR4(7) + grade5*alpha_LR4(11) ...
%               + fem*alpha_LR4(14) + blk*alpha_LR4(18) + hsp*alpha_LR4(19) + SSEInc*alpha_LR4(24) + SSEIns*alpha_LR4(25) ...
%               + NHelpers*alpha_LR4(29) + NoHmIntrnt*alpha_LR4(32) + MoblFrac*alpha_LR4(35) + TabltFrac*alpha_LR4(38) ;  
alphaE_LR4D1  = alphaE_LR4 - D2*alpha_LR4(6) - D3*alpha_LR4(7) ;  
alphaE_LR4D2  = alphaE_LR4 - D2*alpha_LR4(6) - D3*alpha_LR4(7) + alpha_LR4(6) ;  
alphaE_LR4D3  = alphaE_LR4 - D2*alpha_LR4(6) - D3*alpha_LR4(7) + alpha_LR4(7) ;  
subplot(3,1,2)
hold on
[F0,A0] = ecdf(alphaE_LR4D1(D2==0&D3==0)) ;
[F1,A1] = ecdf(alphaE_LR4D1) ;
stairs(-A0,1-F0,'r','LineWidth',1) ;
stairs(-A1,1-F1,'r','LineWidth',4) ; 
[F0,A0] = ecdf(alphaE_LR4D2(D2==1)) ;
[F1,A1] = ecdf(alphaE_LR4D2) ;
stairs(-A0,1-F0,'--b','LineWidth',1) ;
stairs(-A1,1-F1,'--b','LineWidth',4) ; 
[F0,A0] = ecdf(alphaE_LR4D3(D3==1)) ;
[F1,A1] = ecdf(alphaE_LR4D3) ;
stairs(-A0,1-F0,'-.k','LineWidth',1) ;
stairs(-A1,1-F1,'-.k','LineWidth',4) ; 
xlim([0,0.4]) 
xlabel('(-1)\times\alpha_{pi} (productivity production weight)','FontWeight','bold','FontSize',28) 
ylabel('CDFs','FontWeight','bold','FontSize',28) 
legend('District 1 (within, ToT)','District 1 (treatment)','District 2 (within, ToT)','District 2 (treatment)',...
    'District 3 (within, ToT)','District 3 (treatment)','Location','EastOutside','FontWeight','bold','FontSize',18) 
box on
grid on
disp('************************************KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS************************************')
disp('************************************KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS************************************')
disp('KOLMOGOROV-SMIRNOV TESTS FOR DIFFERENCES OF CDFS OF HETEROGENEOUS TREATMENT EFFECTS: SCHOOL ON PRODUCTIVITY')
disp('D1 vs D2 treatments:')
[h,p,ks2stat]=kstest2(alphaE_LR4D1,alphaE_LR4D2) %#ok
disp('D1 vs D3 treatments:')
[h,p,ks2stat]=kstest2(alphaE_LR4D1,alphaE_LR4D3) %#ok
  

% alphaL_LR4  = alpha_LR4(3)*ones(length(logPreScore),1) + D2*alpha_LR4(8) + D3*alpha_LR4(9) + grade5*alpha_LR4(12) ...
%               + fem*alpha_LR4(15) + blk*alpha_LR4(20) + hsp*alpha_LR4(21) + SSEInc*alpha_LR4(26) + SSEIns*alpha_LR4(27) ...
%               + NHelpers*alpha_LR4(30) + NoHmIntrnt*alpha_LR4(33) + MoblFrac*alpha_LR4(36) + TabltFrac*alpha_LR4(39) ; 
alphaL_LR4D1  = alphaL_LR4 - D2*alpha_LR4(8) - D3*alpha_LR4(9) ; 
alphaL_LR4D2  = alphaL_LR4 - D2*alpha_LR4(8) - D3*alpha_LR4(9) + alpha_LR4(8) ; 
alphaL_LR4D3  = alphaL_LR4 - D2*alpha_LR4(8) - D3*alpha_LR4(9) + alpha_LR4(9) ; 
subplot(3,1,3)
hold on
[F0,A0] = ecdf(alphaL_LR4D1(D2==0&D3==0)) ;
[F1,A1] = ecdf(alphaL_LR4D1) ;
stairs(-A0,1-F0,'r','LineWidth',1) ;
stairs(-A1,1-F1,'r','LineWidth',4) ;
[F0,A0] = ecdf(alphaL_LR4D2(D2==1)) ;
[F1,A1] = ecdf(alphaL_LR4D2) ;
stairs(-A0,1-F0,'--b','LineWidth',1) ;
stairs(-A1,1-F1,'--b','LineWidth',4) ;
[F0,A0] = ecdf(alphaL_LR4D3(D3==1)) ;
[F1,A1] = ecdf(alphaL_LR4D3) ;
stairs(-A0,1-F0,'-.k','LineWidth',1) ;
stairs(-A1,1-F1,'-.k','LineWidth',4) ;
xlim([0,0.075])
xlabel('(-1)\times\alpha_{mi} (motivation production weight)','FontWeight','bold','FontSize',28)
ylabel('CDFs','FontWeight','bold','FontSize',28) 
legend('District 1 (within, ToT)','District 1 (treatment)','District 2 (within, ToT)','District 2 (treatment)',...
    'District 3 (within, ToT)','District 3 (treatment)','Location','EastOutside','FontWeight','bold','FontSize',18) 
box on
grid on
% plot([0,0],[0,1],'k') 
disp('KOLMOGOROV-SMIRNOV TESTS FOR DIFFERENCES OF CDFS OF HETEROGENEOUS TREATMENT EFFECTS: SCHOOL ON MOTIVATION')
disp('D1 vs D2 treatments:')
[h,p,ks2stat]=kstest2(alphaL_LR4D1,alphaL_LR4D2) %#ok
disp('D1 vs D3 treatments:')
[h,p,ks2stat]=kstest2(alphaL_LR4D1,alphaL_LR4D3) %#ok

% TFP_LR4     = exp( alpha_LR4(1)*ones(length(logPreScore),1) + D2*alpha_LR4(4) + D3*alpha_LR4(5) + grade5*alpha_LR4(10) ...
%               + fem*alpha_LR4(13) + blk*alpha_LR4(16) + hsp*alpha_LR4(17) + SSEInc*alpha_LR4(22) + SSEIns*alpha_LR4(23) ...
%               + NHelpers*alpha_LR4(28)  + NoHmIntrnt*alpha_LR4(31) + MoblFrac*alpha_LR4(34) + TabltFrac*alpha_LR4(37) ) ; 
TFP_LR4D1     = exp( log(TFP_LR4) - D2*alpha_LR4(4) - D3*alpha_LR4(5) ) ;
TFP_LR4D2     = exp( log(TFP_LR4) - D2*alpha_LR4(4) - D3*alpha_LR4(5) + alpha_LR4(4) ) ;
TFP_LR4D3     = exp( log(TFP_LR4) - D2*alpha_LR4(4) - D3*alpha_LR4(5) + alpha_LR4(5) ) ;
subplot(3,1,1)
hold on
[F0,A0] = ecdf(TFP_LR4D1(D2==0&D3==0)) ;
[F1,A1] = ecdf(TFP_LR4D1) ;
stairs(A0,F0,'r','LineWidth',1) ;
stairs(A1,F1,'r','LineWidth',4) ;
[F0,A0] = ecdf(TFP_LR4D2(D2==1)) ;
[F1,A1] = ecdf(TFP_LR4D2) ;
stairs(A0,F0,'--b','LineWidth',1) ;
stairs(A1,F1,'--b','LineWidth',4) ;
[F0,A0] = ecdf(TFP_LR4D3(D3==1)) ;
[F1,A1] = ecdf(TFP_LR4D3) ;
stairs(A0,F0,'-.k','LineWidth',1) ;
stairs(A1,F1,'-.k','LineWidth',4) ;
xlim([10.5,20.5])
xlabel('A_i (TFP)','FontWeight','bold','FontSize',28)
ylabel('CDFs','FontWeight','bold','FontSize',28) 
legend('District 1 (within, ToT)','District 1 (treatment)','District 2 (within, ToT)','District 2 (treatment)',...
    'District 3 (within, ToT)','District 3 (treatment)','Location','EastOutside','FontWeight','bold','FontSize',18)
box on 
grid on
% plot([0,0],[0,1],'k') 
disp('KOLMOGOROV-SMIRNOV TESTS FOR DIFFERENCES OF CDFS OF HETEROGENEOUS TREATMENT EFFECTS: SCHOOL ON TFP')
disp('D1 vs D2 treatments:')
[h,p,ks2stat]=kstest2(TFP_LR4D1,TFP_LR4D2) %#ok
disp('D1 vs D3 treatments:')
[h,p,ks2stat]=kstest2(TFP_LR4D1,TFP_LR4D3) %#ok
disp('************************************KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS************************************')
disp('************************************KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS-KS************************************')
%%
pre = pretemp ;